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Abstract: We present a dynamical study of the double parton distribution in impact 
parameter space, which enters into the double scattering cross section in hadronic collisions. 
This distribution is analogous to the generalized parton densities in momentum space. We 
use the Lund Dipole Cascade model, presented in earlier articles, which is based on BFKL 
evolution including essential higher order corrections and saturation effects. As result 
we find large correlation effects, which break the factorization of the double scattering 
process. At small transverse separation we see the development of "hot spots", which 
become stronger with increasing Q 2 . At smaller x-values the distribution widens, consistent 
with the shrinking of the diffractive peak in elastic scattering. The dependence on Q 2 
is, however, significantly stronger than the dependence on x, which has implications for 
extrapolations to LHC, e.g. for results for underlying events associated with the production 
of new heavy particles. 
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1. Introduction 

Analyses of high energy pp collisions show that multiple hard parton collisions are quite 
common. Four-jet events, in which the jets balance each other pairwise, were observed 
at y/s = 63 GeV by the AFS collaboration at the CERN ISR [1], and at the Tevatron 
events with four jets, or with three jets + 7, have been observed by the CDF [2, 3] and 
DO [4,5] collaborations. These results also show that the hard subcollisions are correlated, 
and double interactions occur with a larger probability than expected for uncorrelated hard 
interactions. 

Knowing the cross section for double subcollisions is important for an estimate of the 
underlying event. A good understanding of the correlations is therefore also essential for 
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the interpretation of signals for new physics at the LHC. The aim of this paper is to study 
what kind of dynamical correlations in momentum space and impact parameter space are to 
be expected at higher energies, as a result of parton evolution to small x and of saturation 
effects. 

Correlations which follow from energy and flavour conservation have been discussed 
e.g. in refs. [6-8], and from an impact parameter picture in ref. [9]. Nontrivial correlations 
which are consequences of DGLAP evolution have been studied in refs. [7, 10, 11]. In 
ref. [12] it is shown that this implies an increasing correlation for higher Q 2 . However, in 
these references the connection between correlations in momentum and 6-space have not 
been studied. Correlation effects also follow from fluctuations in the number of particles in 
the parton cascades, which is discussed in refs. [6,7]. For related problems in connection 
with production of two heavy gauge bosons, see e.g. refs. [13,14]. 

In the model by Sjostrand and van Zijl [8], later modified in ref. [15] and implemented in 
the PYTHIA event generator [16-18], it is assumed that the dependence of the double parton 
density on the kinematic variables (x and Q 2 ) and the separation in impact parameter space 
(b) factorizes. The correlation is described by a distribution in the separation b, where a 
denser central region implies that central collisions have on average more hard interactions, 
and peripheral collisions fewer. The relative shape of this distribution is kept constant, but 
the width is scaled with energy, to accommodate the growth of the total (non-diffractive) 
and the multiple interaction cross sections with increasing energy. For the dependence on 
x and Q 2 effects of energy and flavour conservation are taken into account. In a recent 
modification of the model, the shape is also varying with x [19], in a way which gives smaller 
correlations for smaller x. At a fixed energy, higher Q 2 is related to larger x, which implies 
that the production of heavy mass objects is more common in high multiplicity events. 
Similar approaches are also implemented in the HERWIG++ event generator [20-22]. 

While it is straight forward to compare the models in PYTHIA and HERWIG with exper- 
imental data on double parton scattering, it is not trivial to disentangle exactly the sources 
of possible correlations in the underlying models. This makes it difficult to extrapolate to 
higher energy and/or higher Q 2 . Here we will instead use a detailed dynamical model for 
parton evolution, for a study of the correlations between gluons inside hadrons. 

In a series of papers we have presented a model based on BFKL [23, 24] evolution and 
saturation, which well reproduces data on total, elastic and diffractive cross sections in 
DIS and pp collisions [25-29]. It is based on Mueller's dipole cascade model [30-32], but 
includes also non-leading effects from energy conservation and running coupling, as well as 
confinement effects and saturation within the evolution. The model is implemented in a 
Monte Carlo (MC) program, which makes it easy to study in detail the types of correlations 
and fluctuations in x and Q 2 as well as in impact parameter space, which are consequences 
of the parton evolution. 

The correlations come from two different sources. First the distribution in transverse 
separation will vary with energy and virtuality, as the impact parameter profile widens at 
higher energy, but at the same time the development of "hot spots" makes the effective 
double scattering profile more narrow. Secondly, BFKL dynamics implies large fluctuations 
in the cascade evolution. This implies that double scattering is more probable in events with 
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a higher than average number of partons. This effect is neglected in most phenomenological 
analyses, which are based on average parton densities. 

The outline of this paper is as follows. First we will go through experimental aspects 
and the theoretical framework for double parton scattering in section ^. In section |3| we will 
then briefly describe the dipole model in the DIPSY Monte Carlo, followed by a description 
in section || of how we will use this model to investigate correlations in the double parton 
densities. The results are presented in section |B| followed by conclusions and outlook in 
section ||. 



2. Double Parton Scattering and Double Parton Distributions 
2.1 Experimental results 

A measure of the correlation is given by the quantity cr e g defined by the relation 

Here a® A B ^ is the cross section for the two hard processes A and B, and the 
corresponding single inclusive cross sections, and (1 + 5ab) _1 is a symmetry factor equal to 
1/2 if A = B. If the hard interactions were uncorrelated, cr e fr would be equal to the total 
non-diffractive cross section. The CDF and DO measurements give instead cr e g ~ 15 mb, 
which thus is significantly smaller. 

As mentioned in the introduction, experimental results on double parton interactions 
have been published from the ISR and the Tevatron. (The UA2 collaboration at the 
CERN SppS collider has presented a lower limit for cr e ff.) The ISR experiment [1] studied 
events with four jets with a summed transverse energy above 29 GeV in pp collisions at 
y/s = 63 GeV. The observed rate corresponds to <j e g = 5 mb, indicating a very strong 
correlation between the partons. As the sum of E± for the four jets is about 30 GeV, 
the typical x-values for the partons is about 0.25. These partons must be dominated by 
valence quarks, and therefore not representative for the partons involved at higher energies 
and smaller x. 

The CDF collaboration studied 4-jet events with J2p±j c t > 140 GeV at y/s = 1.8 TeV 
[2], which implies x- values ~ 0.04. The result obtained was a c s = 12.1^° 4 7 mb. CDF 
also studied events with 7+3jets [3]. This gave a more clear signal than the 4-jet events, 
and the result is <r e ff = 14.5 ± 1.7^ 2 '3 m b- The photons had -Ej_ 7 > 16 GeV, and the jets 
E±j e t > 5 GeV, and thus the x-values are smaller. No dependence on the Feynman scaling 
variable for the two pairs was observed within the ranges 0.01 < x_f(t + jet) < 0.4 and 
0.002 < Xi?(dijet) < 0.2. We note, however, that the 4-jet and 7+3jet processes are not 
equivalent, as the production of a photon are particularly sensitive to quark distributions. 

While the DO 4-jet signal is rather weak [4], this collaboration also measures a clear 
signal for 7+3jets [5]. The transverse momenta, and thus the parton x-values, are larger 
than those in the CDF analysis; p±^ > 60 GeV and p±j e t2 > 15 GeV. The result obtained 
is <7 e g = 16.4 ± 0.3 ± 2.3 mb. DO also studied the dependence on p± of the second jet. Here 
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<7 e ff dropped by 24% when p±_ increased from 17.5 to 27.5 GeV, but this variation was fully 
within the errors of the measurement. 



2.2 Formalism 

Following ref. [11] we define the "double parton distribution" Tij(xi,X2,b; Q\, Q 2 ), which 
describes the inclusive density distribution for finding a parton of type i with energy fraction 
x\ at scale Q 2 , together with a parton of type j with energy fraction X2 at scale Q 2 , an d 
with the two partons separated by a transverse distance b. The distributions are via 
a Fourier transformation related to the "two-parton generalized parton distributions" in 
transverse momentum space, D(x%,X2, Q 2 , Q|' studied by Blok et al. [33]. 

Assuming factorisation of two hard subprocesses A and B, the cross section for double 
scattering is given by 1 



rrf\ ^— ^2 j Fij(xi,X2,bQl,Q2)^(x 1 ,Xi)&f l (x2,x' 2 ) 

xTki(xi,x 2 ,b; Ql, Q\)dxidx2dx' l dx' 2 d 2 b. (2.2) 



fWl + SAB..,,. 

t,3,k,l 



Here a is the cross section for a parton- level subprocess. 

The approximation used in the event generators PYTHIA and HERWIG means that, 
apart from effects of energy and flavour conservation (which are quite important), T fac- 
torizes in the form 

T ij (x 1 ,x 2 ,b;QlQ 2 2 ) ^ D i (x l ,Ql)D\x 2 ,Ql)F{b;s). (2.3) 

Here D l denotes the single parton distribution for parton type i, and F(b; s) describes the 
distribution in the transverse separation between the two partons, and is assumed to be 
independent of X{ and Q 2 , and to be the same for quarks and gluons. (As mentioned in 
the introduction, in a recent option in PYTHIA, F depends also on x [19].) The width is 
adjusted by tuning the model to reproduce the total cross section and multiple interactions, 
and thus depends on the collision energy y/s, as indicated in eq. ( |2.3|) . 

As mentioned in the introduction, the relation in eq. ( |2.3| ) is not consistent with 
DGLAP evolution. Assume that the distribution of two gluons satisfies the factorization 
relation in (2.3) for some specific values of x\, Q\,X2, Q\- Evolving to higher virtualities, 



gluon 1 and gluon 2 can form separate cascades by emitting new gluons. A pair of gluons 
can then be formed, either with one from each of these cascades, or by two gluons from 
the same cascade. The latter contribution will then necessarily break the factorization 
relation [11]. Gaunt and Stirling assume instead the weaker relation 

T ij (x 1 ,X2,b;QlQ 2 2 ) = DV(x 1 ,x 2 ;QlQ 2 2)Fi(b). (2.4) 

The relations in eq. ( |2.2| ) and eq. (|2.3D or ( |2.4| ) imply that the effective cross section in 
Q2.1| ) is determined by the relation 

-l 



CcfT 



/ d 2 b(F(b)f 



(2.5) 



1 This result is true if the parton-parton scattering is local in impact parameter space, which ought to 
be the case for hard collisions with large Q 2 -values. 



- 4 - 



In the following we want to use the Lund Dipole Cascade model to study the corre- 
lations and fluctuations, which follow from the parton evolution in a proton, and see how 
well the approximations in eqs. ( |2.3[ ) and ( |2.4| ) are satisfied. We note, however, that as our 
model is based on the BFKL evolution, it contains only gluons and should only be trusted 
at small x-values where the gluons dominate. 

We define the distribution F(b; x\, 22, Q 2 , Q2) by the relation 



T(x 1 ,x 2 ,b; Qx,Ql) 



D(x 1 ,Q 2 1 )D(x 2 , Ql)F{b- X1 ,X2, Ql Ql). 



(2.6) 



Thus F is a density in transverse coordinate space b, which may depend on all four variables 
x\, X2, Qi, and Q 2 , and which contains all information about correlations between the two 
partons. In case e.g. some kind of "hot spots" develop for small x and/or large Q 2 , this 
will show up as an increase in F for small b- values. With this definition F is related to the 
double scattering cross section via eq. (|2.1|) and the relation 



Ccff 



with the constraint 



d 2 bF{b- Xl ,x 2 , Ql Ql)F(b; x[,x 2 , QlQl) 



Ql 

x\x\ 



Ql 

x 2 x' 2 ' 



For hard subcollisions at midrapidity we have x\ ~ x\ and X2 
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(2.7) 



(2.8) 



and thus recover 



eq. ( |2.5| ), with the difference that cr e ff may now depend on the variables xi and Qf (with 
Qi/ X i = Q\l x \ = s )- We n °te, however, that for subcollisions away from midrapidity we 
get different x-values in the two F-distributions in eq. ( p.7[ ). This feature will be further 
discussed in sec. ||[ 



2.3 Correlations 

There are two different sources for correlations between the partons, which both give con- 
tributions to a e fi. 



2.3.1 Distribution in impact parameter space 

More central collisions are expected to have more hard subcollisions than peripheral col- 
lisions. This feature causes correlations, which depend on the matter distribution within 
the colliding protons. To understand the effect of the matter distribution, we here study a 
simplified version of the model implemented in the PYTHIA event generator. A discussion 
of the results in our model, and a comparison with more refined versions of PYTHIA, are 
presented in sec. ||. 

Assume that the parton density inside a proton is given by a Gaussian distribution 
p oc exp(— r 2 /a 2 ). Integrating over the longitudinal coordinate z gives a factor \fn a inde- 
pendent of the impact parameter b = \J x 1 + y 2 . Thus the density in impact parameter 
space is also given by 

poc exp(-6 2 /a 2 )- (2-9) 
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The two parton density is then given by 

F(b) oc J d 2 rp(r)p(r - b) oc exp(-6 2 /2a 2 ). (2.10) 

The F-distribution is normalized to 1, and has therefore a normalization constant equal to 
(-7r2a 2 ) -1 . This implies that 

-i 



Ceff 



d 2 bF 2 ) =i-Ka 2 . (2.11) 



In PYTHIA it is assumed that for two colliding protons, the average number of hard 
subcollisions is proportional to the overlap of the two distributions. For a collision at 
impact parameter 6, this overlap is given by 

0(b) = J d 2 rp(r)p(r - b) oc exp(-6 2 /2a 2 ). (2.12) 



We note that the integral in eq. ( 2.12j ) is exactly the same as the one determining F in 



eq. ( |2,10| ). The distribution F is normalized to 1, and we therefore find 

F(b) = 0(b)/ J d 2 bO(b). (2.13) 

Since in PYTHIA the average number of hard subcollisions at a fixed impact parameter, 
n(b), is proportional to the overlap, 0(b), it can be written as 

n(b) oc 0(b) ■ a, (2.14) 

where the cross section for parton-parton collisions is approximated by &S^(h). The 
number of subcollisions is assumed to be given by a Poisson distribution with average nib), 
and a non-diffr active event is obtained if n ^ 0. This implies that the total non-diffractive 
cross section is given by 

ctnd = / d 2 b(l-e- n ^), (2.15) 



and the average number of subcollisions is 

, v fd 2 bn(b) . 

(n) = V -L. (2.16) 

In PYTHIA a depends on a parameter p±o, which acts as a smooth cutoff for small 
p± in parton-parton scattering. We have therefore two adjustable parameters, a and pj_o, 
which can be used to fit cxnd and (n). A result from such a fit is presented in ref. [19], with 
the result a ~ 0.48 fm, giving <7 e ff ~ 28 mb, at the Tevatron and cr e ff ~ 34 mb at LHC. 

These values could be regarded as a kind of benchmarks. If the evolution gives smaller 
regions with higher parton density, the width of the F-distribution will be wider, and <7 e ff 
will be smaller. We here note that, although the single Gaussian fit gives the proper number 
of subcollisions, the correlations appear to be too week. More successful tunes to data have 
a distribution p (and thus also O), which is given by a sum of two Gaussian. This enhances 
small and large b- values, and gives a somewhat smaller cr e fj. In a recent modification to 
the model [19], the width of the density p is assumed to vary with x, which also has the 
effect of enhancing small and large b- values increasing the correlations and reducing cr e fj. 



-6- 



/ 



d %F{b) = y^rzA - ( 2 - 18 ) 



2.3.2 Fluctuations in the parton cascade 

Another source for correlations is coming from fluctuations within the cascades, and affects 
a c g even if the double parton distribution factorizes as in eq. (|2.3| ) . This effect is discussed 
in refs. [6] and [7], but is normally neglected in phenomenological analyses. Double hard 
scattering is more likely in collisions with protons that have more than the average number 
of partons. A measure of this effect is given by the integral of F as follows: 

The cascade evolution is a random process, which can lead to different partonic states. 
We label the states by the parameter n, and the parton distribution in state n is denoted 
D n (x, Q 2 ). The probability to obtain this state is denoted P n , with ^ P n = 1. For a fixed 
state n the density of partons in a small interval 5xi around Xi at resolution Qf, is given 
by D n (xi, Qf)5xi. As we have a single state, the existence of a parton at x± is fixed, and it 
is not possible to define a correlation with the existence of another parton at x%. Therefore 
Dn{xi,Ql,X2, Q2) = D n (xi, Q\) ■ D n (x2, Q 2 )- Averaging over all states n we thus get 

d 2 bF(b) _ EnPnDn(xi,Ql)D n (x 2 ,Q 2 2 ) (D(x U Q\)D{x 2 , Qj)) 

EnPnDn(xi,Ql)-^: n PnD n (x2,Q 2 2 ) {D{ X% , Q 2 )) {D(x 2 , Q 2 )) ' K ' ' 

Here (...) denotes an average over the different cascades n. For the special case of Q 2 = Q 2 
and x\ = x 2 , we get 

(D 2 (x,Q 2 )} 
{D(x,Q 2 ))< 

Thus the integral of F equals 1 + V/(D(x, Q 2 )} 2 , with the variance V equal to the square 
of the width of the distribution. We see that the fluctuations imply that the F-distribution 
is enhanced, which also implies a larger correlation and a smaller cr e g. When presenting 
our results in section ||, we will also include values for the integral of F, in order to 
facilitate the interpretation of the results. Note that in the Gaunt-Stirling approximation 
in eq. (|2,4|), F is normalized to 1, and all correlation effects are included in the function 
D^(x 1 ,x 2 ;Q 2 ,Q 2 ). 

3. The Lund Dipole Cascade Model 
3.1 Mueller's dipole cascade 

Mueller's dipole cascade model [30-32] is a formulation of the leading logarithmic (LL) 
BFKL evolution in transverse coordinate space. Gluon radiation from the colour charge 
in a parent quark or gluon is screened by the accompanying anti-charge in the colour 
dipole. This suppresses emissions at large transverse separation, which corresponds to the 
suppression of small k± in BFKL. For a dipole (x, y) the probability per unit rapidity (Y) 
for emission of a gluon at transverse position z is given by 

dV a o fx — y) 2 _ 3a s 

— = —d z- rj, with a = . (3.1) 

dY 2tt (x — z) z (z — y) z it 

This emission implies that the dipole is split into two dipoles, which (in the large N c limit) 
emit new gluons independently. The result is a cascade, where the number of dipoles grows 
exponentially with Y. 
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In a high energy collision, the dipole cascades in the projectile and the target are 
evolved from their rest frames to the rapidities they will have in the specific Lorentz frame 
chosen for the analysis. The scattering probability between two elementary colour dipoles 
with coordinates and (xj,yj) in the projectile and the target respectively, is given 

by 2/jj, where in the Born approximation 

{ Xi - Xj )2 {y .-y.)2)\ ■ ^ > 

The optical theorem then implies that the elastic amplitude for dipole i scattering off dipole 
j is given by fy. Summing over i and j gives the one-pomeron elastic amplitude 

F = Y,fir (3-3) 

The growth in the number of dipoles also implies a strong growth for the interaction 
probability, but the total scattering probability is kept below 1 by the possibility to have 
multiple dipole-dipole subcollisions in a single event. In the eikonal approximation the 
unitarized elastic amplitude is given by the exponentiated expression 

T(b) = 1 - e~ F , (3.4) 

and the total, elastic, and diffr active cross sections are given by 

da tot /d 2 b = (2T), 
da cl /d 2 b = (T) 2 , 
da diSex /d 2 b = (T 2 ) - (T) 2 . (3.5) 

3.2 The Lund dipole cascade model 

In refs. [25,26,28] we describe a modification of Mueller's cascade model with the following 
features: 

• It includes essential NLL BFKL effects. 

• It includes non-linear effects in the evolution. 

• It includes effects of confinement. 

The model also includes a simple model for the proton wavefunction, and is imple- 
mented in a Monte Carlo simulation program called DIPSY. As discussed in the cited refer- 
ences, the model is able to describe a wide range of observables in DIS and pp scattering, 
with very few parameters. 

3.2.1 NLL effects 

The NLL corrections to BFKL evolution have three major sources [34]: 
The running coupling: 

This is relatively easily included in a MC simulation process. The scale in the running 
coupling is chosen as the largest transverse momentum in the vertex [35]. 

Non-singular terms in the splitting function: 
These terms suppress large z-values in the individual parton branchings, and prevent the 



fij — f{ x 'hyi\ x jiyj) — o 



log 
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daughter from being faster than her recoiling parent. Most of this effect is taken care of by 
including energy-momentum conservation in the evolution. This is effectively taken into 
account by associating a dipole of transverse size r with a transverse momentum k± = 1/r, 
and demanding conservation of the light-cone momentum p + in every step in the evolution. 
This gives an effective cutoff for small dipoles, which also eliminates the numerical problems 
encountered in the MC implementation by Mueller and Salam [36]. 

Projectile-target symmetry: 
This is also called energy scale terms, and is essentially equivalent to the so called con- 
sistency constraint [37]. This effect is taken into account by conservation of both positive 
and negative light-cone momentum components, p+ and p_. The treatment of these effects 
includes also effects beyond NLL, in a way similar to the treatment by Salam in ref. [34]. 
Therefore the power A e g, determining the growth for small x, does not turn negative for 
large values of a s . 

3.2.2 Non-linear effects and saturation 

As mentioned above, dipole loops (or equivalently pomeron loops) are not included in 
Mueller's cascade model, if they occur within the evolution. They are only included if they 
are cut in the Lorentz frame used in the calculations, as a result of multiple scattering in this 
frame. The result is therefore not frame independent. (The situation is similar in the Color 
Glass Condensate [38-40] or the JIMWLK [41,42] equations.) As for dipole scattering the 
probability for such loops is given by a s , and therefore formally colour suppressed compared 
to dipole splitting, which is proportional to a = N c a s /TT. These loops are therefore related 
to the probability that two dipoles have the same colour. Two dipoles with the same 
colour form a quadrupole field. Such a field may be better approximated by two dipoles 
formed by the closest colour-anticolour charges. This corresponds to a recoupling of the 
colour dipole chains, favouring the formation of small dipoles. We call this process a dipole 
"swing". The swing gives rise to loops within the cascades, and makes the cross section 
frame independent to a good approximation. We note that a similar effect would also be 
obtained from gluon exchange between the two dipoles. 

In the MC implementation each dipole is assigned one of colour indices, and dipoles 
with the same colour index are allowed to recouple [26]. The weight for the recoupling is 
assumed to be proportional to (rfr^/ir^rj), where r\ and r2 are the sizes of the original 
dipoles and and are the sizes of the recoupled dipoles. Dipoles with the same colour 
are allowed to swing back and forth, which results in an equilibrium, where the smaller 
dipoles have a larger weight. We note that in this formulation the number of dipoles is 
not reduced, and the saturation effect is obtained because the smaller dipoles have smaller 
cross sections. Thus in an evolution in momentum space the swing would not correspond to 
an absorption of gluons below the saturation line k\ = Q 2 s {x); it would rather correspond 
to lifting the gluons to higher k± above this line. Although this mechanism does not 
give an explicitly frame independent result, MC simulations show that it is a very good 
approximation. 
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3.2.3 Confinement effects 



Confinement effects are included via an effective gluon mass, which gives an exponential 
suppression for very large dipoles [27]. This prevents the proton from growing too fast in 
transverse size, and is also essential to satisfy Froisart's bound at high energies [43]. 

3.2.4 Initial dipole configurations 

In DIS an initial photon is split into a qq pair, and for larger Q 2 the wavefunction for 
a virtual photon can be determined perturbatively. The internal structure of the proton 
is, however, governed by soft QCD, and is not possible to calculate perturbatively. In 
our model it is represented by an equilateral triangle formed by three dipoles, and with a 
radius of 3 GeV -1 f» 0.6 fm. The model should be used at low x, and when the system is 
evolved over a large rapidity range the observable results depend only weakly on the exact 
configuration of the initial dipoles, or whether the charges are treated as (anti)quarks or 
gluons. 

4. Application to double parton distributions 

In principle we could estimate the gluon density in the proton from the cross section 
for 7*p collisions. The photon would then be treated as a superposition of dipoles of 
varying sizes, with weights determined by QED. However, in order to more easily isolate the 
correlations and fluctuations in the proton, without the complications from the additional 
fluctuations in the photon wave functions, we prefer instead to calculate the cross section 
for a dipole with a fixed transverse size. The photon coupling to a qq pair contains a 
factor (z(l - z)) 2 ~ i K 2 {^/z{\ - z) Qr). Here K, L are generalized Bessel functions, r is the 
transverse separation between the quark and the antiquark, z is the fractional energy 
taken by the quark, and the index i is 1 for transverse and for longitudinal photons. 
Important contributions are obtained for z ~ 0.5, and since the Bessel functions fall off 
exponentially when the argument is larger than 1, characteristic r-values are given by the 
relation r « 2/Q. 

To determine the parton distributions D(x,Q 2 ) we thus calculate the cross section for 
scattering of a single dipole of size r = 2/Q colliding with a proton. The projectile dipole 
has a large cross section when colliding with equally large or larger dipoles in the proton, 
while the interaction with smaller dipoles is suppressed. The dipoles in the proton are also 
connected to gluons with k± of the order 1/r, and therefore the D-distributions obtained 
correspond to the integrated gluon distributions. 

The double parton density r(xi, X2, b; Q 2 , Q 2 .) is in the same way proportional to the 
cross section, a^' 2 \ for simultaneous scattering of two dipoles with size V{ = 2/Qi, sepa- 
rated by a distance b, and colliding with a proton. Thus also T is defined as an integrated 
density. From these relations we conclude that the F-distribution is directly related to the 
ratio between the corresponding cross sections, and in an obvious notation we have 



F(b; Xl ,X2,Ql,Ql) = 



r(i> 2 ) 



a 



(1-2) 



(4.1) 



dQ-)dw 



a 



(1)<t(2)- 



- 10 - 




Figure 1: Correlation function F(b) corresponding to central subcollisions at t/s ~ 1.5 TeV, and 
three combinations of Q\ and Q\. The x-values corresponding to Q 2 = 10 and 1000 GeV 2 are 10~ 3 
and 10 -2 respectively. 



In principle the double scattering cross section, and thus also F, depends in the dipole 
model on the three vectors, the dipole sizes ri, r2, and the distance between them b. This 
implies that the estimate of <7 e ff should be given by a 6-dimensional integral over these 
variables. For the results presented in the next section we have calculated F keeping the 
separation between two gluons constant equal to b, but averaging over the directions of the 
dipoles, ri and r2. To check this approximation we have also calculated the result obtained 
when keeping the distance between the centers of the dipole fixed. The results presented 
in sec. [5] show that the correlations obtained are insensitive to how the averages are taken. 

5. Results 

In this section we show some results relevant for pp collisions at ^/s ~ 1.5 and 15 TeV, 
qualitatively representing Tevatron and LHC energies. We calculate the cross sections for 
one or two dipoles against a proton, where the rapidity separation between the projectile 
and the target is given by y = ln(l/xj) + ln(Qj/2). Thus for subcollisions at midrapidity, 
with y = ln-^/i, we have x = Q/(2^/s). 

5.1 Subcollisions at midrapidity 

In figs. [I] and || we show the Q 2 -dependence for x-values corresponding to central production 
at a fixed energy. Fig. [l] shows the result for ^fs « 1.5 TeV and three combinations 2 of 
Qi and Q\. The corresponding x-values are given by x = Q/(2^). Fig. || shows similar 
results for yfs « 15 TeV. 

We see that the distribution is not well described by a Gaussian. Instead there is 
an almost exponential dependence in the range 0.2 < b < 0.6 fm. For larger b- values the 

2 The chosen values of Q 2 of 10, 10 3 and 10 5 GeV 2 , correspond to dipole sizes of w 0.13, 0.013 and 
0.0013 fm respectively. 
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Figure 2: Correlation function F(b) corresponding to central subcollisions at ^fs ~ 15 TeV. The 
curves correspond to Q\ = Q\ = 10, 10 3 , andl0 5 GeV 2 . The corresponding z-values are 10 -4 , 
10" 3 , and 1CT 2 . 



distribution drops faster , and this effect is stronger for higher Q 2 . The position of the 
break is related to the size of the initial proton configuration in the model. For small 
b- values a spike is developing, growing stronger with increasing Q 2 . This is associated with 
a faster drop for larger b for high Q 2 . Results with one softer and one harder subcollision 
lie (as expected) in between the results for two soft or two hard collisions, but closer to the 
result for two softer subcollisions (not shown for yfs ~ 15 TeV). 

In fig. ^ we show the energy dependence by comparing the results for Q 2 = 10 3 GeV 2 
at the two different energies. We see here that the peaks at small b are almost identical, 
but at the higher energy the distribution becomes a little wider, with a slightly larger tail 
out to large 6- values. 

The corresponding values for <7 e fj = (J d 2 b F(b) 2 )^ 1 and J d 2 bF(b) were calculated 
numerically, and are presented in table |]. The rather slow fall off for F at large b gives 
a large a e g for small Q 2 , but the more narrow distributions for higher Q 2 imply a very 
strong Q 2 -dependence. At 15 TeV a e g drops by a factor 2 when Q 2 is changed from 10 to 
10 2 GeV 2 . For fixed Q 2 the wider distribution at higher energy, implies a slightly larger 
cr c fj. The variation with y^s is, however, much weaker than the variation with Q 2 . 

We also note that the effect of fluctuations in the cascades, represented by the difference 
from 1 of J d 2 bF, is of the order 5-10% (contributing 10-20% to cr c g), and is smaller for 
large Q 2 . 

5.2 Sub-collisions off midrapidity 

As mentioned in sec. |2.2| , for subcollisions away from midrapidity, the effective cross section 
is given by eq. ( |2.7D , which contains a product of two ^-distributions with different x- values. 
In fig. ID we show results relevant for two subcollisions with Q 2 = 10 GeV 2 at 1.5 TeV. One 
collision is at rapidity zero, with x\ = x' x = 0.001, while the other subcollision has X2 = 0.01 
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Q^=Q2=10 3 GeV 2 \ls=1.5TeV 




1Q-4 I I I I I I I I I I I I I I I I I I I I I 

0.5 1 1.5 2 



b[fm] 

Figure 3: The energy dependence of the correlation functions F(b). The curves show F for 
Ql = Ql = 10 3 GeV 2 , solid line for y/s « 1.5 and dashed line for y/s « 15 TeV. 



QlQl [GeV 2 ],x 1 ,x 2 


a eS [mb] 


IF 


1.5 TeV, midrapidity 






10 10 0.001 0.001 


35.3 


1.09 


10 10 3 0.001 0.01 


31.0 


1.07 


10 3 10 3 0.01 0.01 


23.1 


1.06 


15 TeV, midrapidity 






10 10 0.0001 0.0001 


40.4 


1.11 


10 3 10 3 0.001 0.001 


26.3 


1.07 


10 3 10 5 0.001 0.01 


24.2 


1.05 


10 5 10 5 0.01 0.01 


19.6 


1.03 


1.5 TeV, y pair2 = 2.3 






10 10 0.001 0.0001 
10 10 0.001 0.01 


} 37.1 


1.08 



Table 1: Summary of results for a e s and corresponding integrals of the double distribution func- 
tions. (The numerical uncertainties are about 1%.) 

and x' 2 = 0.0001, which corresponds to a rapidity for the pair of jets (or a produced massive 
particle) given by y pa i r = ]n(x2/x 2 )/2 « 2.3. The corresponding F-distributions are called 
Fs (for small X2 = 0.0001) and Fl (for large X2 = 0.01). We see that in the tail Fs and Fl lie 
on opposite sides of the distribution -F ce ntrai> which is the relevant one for two subcollisions 
at midrapidity (both with x% = x\ = 0.001). As a e g is given by (j d 2 bFsFi)^ 1 , this implies 
that although the ^-distributions vary with x, <7 e ff is approximately independent of y pa ir- 
This result is further illustrated in fig. || which shows the ratio s/F^Fl /-^central- We see 
that this ratio is close to 1 except for small b- values, which have a low weight in the integral 
over d 2 b. As a consequence also the result for <r e g varies very little with rapidity, in this 
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Q;=Q2=10 GeV 2 x,=0.001 x 2 =0.0001 
Q 2 =Q 2 =10GeV 2 x,=0.001 x 2 =0.001 
Q 2 =Q 2 =10GeV 2 x,=0.001 x 2 =0.01 




0.5 



1.5 



b[fm] 



Figure 4: Correlation function F(b) relevant for subcollisions at y pa ir2 = 2.3 for ,/s s» 1.5 TeV and 
Ql = Ql = 10 GeV 2 . F s = F(b; 0.001, 0.0001, 10, 10) and F L = F{b; 0.001, 0.01, 10, 10). (Q 2 in 
GeV 2 ) The distribution F(b; 0.001, 0.001, 10, 10) for two subcollisions at midrapidity is included 
for comparison. 




Figure 5: The ratio VWFs / F ccntrlih Q\ = Q\ = 10 GeV 2 , and ^/s 1.5 TeV. 



case from 35.3 to 37.1 mb when ypair is changed from to 2.3 (see table [l]). Actually, 
for small shifts 5y pa ir, the difference between the product F^Fs and F^ entral must be only 
second order in <% pa ir- This must then also be the case for the corresponding values for 

Ceff- 



5.3 Comparison with experiment 

The result for <r e g is larger than the effective cross section observed for 7+3jet events at the 
Tevatron. We note, however, that these data depend on quark distributions, and not only 
on the gluon distributions. It would be interesting to try to estimate the difference, but this 
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goes beyond the scope of this paper. For Q 2 = 10 3 GeV 2 at 1.5 TeV we find cr e g = 23 mb. 
Within the errors this is in agreement with the 4-jet results from CDF, which however, 
have large uncertainties. It would also be interesting to have, for comparison, the values 
for <7 e ff in the tuned versions of PYTHIA, not only for the less successful single Gaussian 
approximation used in section 2.3.1. 

We also note that there are two observed features, which are well reproduced by our 
model: 

i. The DO collaboration [5] has observed a reduction in <T e fr for growing pj_ of the jet 
pair in 7 + 3 jet-events. Although the errors are large, the central values for <r e ff drop from 
18.2 to 13.9 mb, when ) increases from 17.6 to 27.3 GeV. This agrees qualitatively 
with the variation with Q 2 found in our model. 

ii. The CDF collaboration [3] observes in the same reaction, that a e s is very insensitive 
to a change in the rapidity y pa i r (or equivalently x V p U ) for one of the subcollisions. As 
discussed in sec. |5.2| this is also consistent with our model. We note, however, that this 
result is mainly a consequence of the fact that the product Fl F$ is close to i^entral , and 
does not neccesarily show that F does not vary with x. 

These qualitative agreements with experimental data may indicate that, also if our 
model possibly underestimates the correlations, we may believe that the qualitative fea- 
tures, and variations with Q 2 and y/s, are correct. The indicated strong dependence on 
Q 2 should then be very important for the interpretation of multiple hard interactions at 
LHC. We here note a possible qualitative difference between our model and the model with 
x-dependent matter densities in ref. [19]. Our result shows a larger variation with Q 2 for 
fixed energy, than with x (or energy) for fixed Q 2 . At fixed energy a variation with Q 2 
is equivalent to a variation with x = Q/y/s. However, the two models may give different 
results when data at different energies are compared. We would like to study this question 
more in the future. 



5.4 Comment on the definition of b 

As mentioned in sec. |||, in the dipole model F depends in principle on three vectors: the 
size of the two dipoles (ri, and the separation between them. For the results presented 
above we have defined the separation as the distance between two gluons, and averaged 
over ri and T2. When calculating J d 2 bF(b) this gives exactly the correct result. This 
is, however, not the case for the integral J d 2 bF 2 (b), which determines a e g. When the 
separation is of the same order of magnitude as the size of the dipoles, the averaging does 
not give the correct result. To estimate this error we have changed the definition of b, 
to the separation between the centers of the dipoles, keeping this fixed when averaging 
over angles for ri and T2- The ratio between the two different F-distributions is shown 
in fig. ||. In this example we have Q 2 = Q\ = 10 GeV 2 and yfs ~ 15 TeV. We can here 
see a difference when b is of the same order of magnitude as the dipoles, around 0.1 fm. 
The contribution from small 6-values is, however, suppressed in the integral over d 2 b, and 
as mentioned above J d 2 bF 2 (b) obtains its major contribution for b in the range 0.5-1 fm. 
The difference between the two estimates of <r e g is only about 2.5%. For larger Q 2 the 
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Figure 6: The ratio between the ^-distributions obtained by two different definitions of b. Q\ — 
Ql = 10 GeV 2 , and y/1 w 15 TeV. 

difference is even smaller, as the region where the two 6-definitions give different results 
moves towards smaller 6-values. This effect can therefore be safely neglected. 

6. Conclusions and outlook 

The Lund Dipole Cascade model offers unique possibilities to study the evolution of gluons 
inside hadrons at small x. The formalism is based on BFKL evolution including essential 
higher-order corrections and saturation effects. By following the evolution, emission-by- 
emission, in rapidity and in transverse position, we can investigate the correlations and 
fluctuations of the gluon distribution in great detail. 

We have here concentrated on the double parton distribution, which enters into the 
multiple parton scattering cross section in proton-proton collisions. The double parton 
distribution in transverse coordinate space is the analogy of the generalized parton density 
in momentum space. The correlations can be described by a distribution in impact pa- 
rameter space, F(b; x\,x%,Q\, Q 2 .). In many applications this distribution is assumed to be 
independent of the energy fractions and scales of the partons, and the same for quarks and 
gluons. In the PYTHIA event generator the width of the distribution does vary with energy, 
and in a recent study it is assumed to vary with x. Earlier analyses have demonstrated 
that DGLAP evolution implies non-trivial correlations depending on x and Q 2 , but this 
analysis does not give information about the 6-dependence. In our analysis we find that 
the two-parton correlation function F(b) depends in a non-trivial way on all the kinematic 
variables x%, x%, Q 2 , Q 2 . 

For subcollisions at midrapidity, F is directly connected to the "effective cross section" 
used in experimental analyses, via the relation J d 2 bF 2 = . The result that F depends 
on Xi and Qf implies, however, that for subcollisions away from y = 0, cr fj is determined by 
an integral of a product of two F-functions with different x-values, connected by the relation 
Xix\ = Q 2 /s. Therefore a e s is very insensitive to the rapidity of a hard subcollision; when 
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one of the two F-functions in the product is above, the other is below the .F-distribution 
relevant for midrapidity subcollisions. An experimental result showing a weak dependence 
on rapidity is therefore not a proof that F does not depend on the scaling parameters Xi . 

In our formalism we can also study event-by-event fluctuations in the density of par- 
tons, which are usually neglected in other analyses. Neglecting the fluctuations F satisfies 
the normalization condition J d 2 bF = 1, but when taking them into account we obtain a 
larger value for this integral, and therefore larger two-parton correlations. In our analysis 
this effect contributes 10-20% to the value of a e g. 

By studying the correlation function F, we can see explicitly how the emission of 
high-p^ gluons result in so-called hot spots developing for small b- values at large Q 2 . For 
fixed Q 2 the emission of larger dipoles, related to low-p± gluons, implies that the distribu- 
tion widens as we go down in x. (This effect is related to the shrinking of the diffractive 
peak at higher energies.) As a result a e s is reduced for higher Q 2 at fixed energy, but 
increases at higher energy for fixed Q 2 . We see, however, that the result appears to depend 
much stronger on Q 2 than on x or collision energy. 

It is not straight-forward to compare our results for the effective cross section with data 
from the Tevatron and ISR, which correspond to lower values for a c s and thus stronger 
correlations. The ISR data are obtained for very high x- values, and the more exact Tevatron 
results are for 7 + 3 jet events, which necessarily involves quarks. Our calculations include 
only gluons, and should therefore only be applied at small x. It is also difficult to estimate 
a possible difference between quark and gluon distributions. Our results are, however, 
compatible within errors with the four-jet measurement at CDF, where gluons should be 
dominating, and they show qualitative agreements with the observed rapidity independence 
at CDF and the hints of strong Q 2 -dependence at DO. We also note that our results appear 
to be in qualitative agreement with results from PYTHIA tuned to multiple interactions and 
underlying events. 

Even if it is possible that we overestimate the result for a e s to some degree, we are more 
confident about the variation with x and Q 2 , which should have relevance for extrapolations 
to the LHC. Our model here predicts a very strong decrease of the effective cross section 
with increasing Q 2 (especially when both Q 2 and Q\ increases), while the dependence on 
the total energy (~ 1/x) for fixed Q 2 is rather weak. Also the dependence on the rapidity 
of the two subcollisions is quite weak, which follows as a e g here is given by a product of 
two .F-distribution, one which is smaller and one which is larger than the F-distribution 
relevant for midrapidity subcollisions. 

When comparing to the recent option in PYTHIA, with an x-dependent impact param- 
eter profile [19], our result may be quite similar if we study the dependence on Q 2 at a fixed 
cms energy y/s. In our model the correlations increase as a direct consequence of higher 
Q 2 , and in the PYTHIA model they increase because x is higher for large Q 2 at fixed energy. 
The difference between the two models might therefore show up in the extrapolation of 
Tevatron data to LHC energies. 

It should be noted, however, that the extraction of the effective cross sections from data 
is very hard, and it will still be difficult to compare future LHC results to the numbers 
we have presented here. For this reason it is important to compare to event generator 
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predictions for the experimental observables. Fortunately we have now developed the 
DIPSY MC to also generate exclusive final states [44], which means that in the future we 
can in detail simulate exclusive final state observables related to double parton scattering. 
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